import numpy as np
import matplotlib.pyplot as plt
c1="#009ad6"
c2='#f15a22'
c3='#7fb80e'
c4='#aa2116'
c5='#ffe600'

rho=np.logspace(3,0,61)
tin=np.logspace(0.46568,3.4652,61)
x=np.load('de445j2.npy')
y=np.load('445extra.npy')
tau=x[:,:,0:4]/tin.reshape(61,1,1)
tau2=y[:,:,0:4]/tin.reshape(61,1,1)
sigma=x[:,:,4:8]/0.14

plt.figure(dpi=150)
plt.plot(rho,np.ones(len(rho)),'k')

#plt.fill_between(rho,np.percentile(tau,16,axis=1)[:,3],np.percentile(tau,84,axis=1)[:,3],
#                 color=c3,alpha=0.3,edgecolor="none")


#plt.plot(rho,np.mean(tau,axis=1)[:,2],'--',color=c1,lw=5,label=r'$\tau_{\rm PE}$')
plt.plot(rho,np.mean(tau,axis=1)[:,3],color=c3,lw=5,
         label=r'$\tau_{\rm PM}$ (K17PMm)')
#plt.plot(rho,np.mean(tau,axis=1)[:,1],'-.',color=c4,lw=3,label=r'$\tau_{\rm MAP}$')
#plt.plot(rho,np.mean(tau,axis=1)[:,0],':',color=c2,lw=3,label=r'$\tau_{\rm MLE}$')

#plt.plot(rho,np.mean(tau2,axis=1)[:,2],'--',color='k',lw=3,alpha=0.5,
#         label=r'$\tau_{\rm PE}(\sigma_*)$')
plt.plot(rho,np.mean(tau2,axis=1)[:,3],'k-.',alpha=0.5,
         label=r'$\tau_{\rm PM}$ (K17PMm,$\sigma_*$)')
#plt.plot(rho,np.mean(tau2,axis=1)[:,1],'-.',color='k',alpha=0.5,
#         label=r'$\tau_{\rm MAP}(\sigma_*)$')
#plt.plot(rho,np.mean(tau2,axis=1)[:,0],':',color='k',alpha=0.5,
#         label=r'$\tau_{\rm MLE}(\sigma_*)$')

plt.xscale('log')
plt.xlabel(r'$T/\tau_{\rm in}$',fontsize=20)
plt.ylabel(r'$\tau_{\rm out}^{\rm mean}/\tau_{\rm in}$',fontsize=20)
plt.text(2,0.2,'mean',fontsize=20)
#plt.title(r'K17 prior (N = 60, $\sigma_e=\sigma_{\rm SDSS}$)',fontsize=20)
plt.title(r'K17 prior (N = 445, $\sigma_{\rm e}=\sigma_{\rm OGLE}$)',fontsize=20)
plt.xlim(1,1000)
plt.grid()
plt.ylim(0.1,1.5)
#plt.legend(fontsize=13,loc=4,ncol=2)
plt.legend(fontsize=15,loc=4)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.show()
'''
plt.figure(dpi=150)
plt.plot(rho,np.ones(len(rho)),'k')

plt.fill_between(rho,np.percentile(tau,16,axis=1)[:,2],np.percentile(tau,84,axis=1)[:,2],
                 color=c1,alpha=0.3,edgecolor="none")

plt.plot(rho,np.median(tau,axis=1)[:,2],'--',color=c1,lw=5,label=r'$\tau_{\rm PE}$ (S21PE50)')
plt.plot(rho,np.median(tau,axis=1)[:,3],color=c3,lw=5,label=r'$\tau_{\rm PM}$')
plt.plot(rho,np.median(tau,axis=1)[:,1],'-.',color=c4,lw=3,label=r'$\tau_{\rm MAP}$')
plt.plot(rho,np.median(tau,axis=1)[:,0],':',color=c2,lw=3,label=r'$\tau_{\rm MLE}$')

plt.plot(rho,np.median(tau2,axis=1)[:,2],'--',color='k',lw=3,alpha=0.5,
         label=r'$\tau_{\rm PE}(\sigma_*)$')
plt.plot(rho,np.median(tau2,axis=1)[:,3],color='k',alpha=0.5,
         label=r'$\tau_{\rm PM}(\sigma_*)$')
plt.plot(rho,np.median(tau2,axis=1)[:,1],'-.',color='k',alpha=0.5,
         label=r'$\tau_{\rm MAP}(\sigma_*)$')
plt.plot(rho,np.median(tau2,axis=1)[:,0],':',color='k',alpha=0.5,
         label=r'$\tau_{\rm MLE}(\sigma_*)$')

plt.xscale('log')
plt.xlabel(r'$T/\tau_{\rm in}$',fontsize=20)
plt.ylabel(r'$\tau_{\rm out}^{\rm median}/\tau_{\rm in}$',fontsize=20)
plt.text(2,0.2,'median',fontsize=20)
#plt.title(r'K17 prior (N = 60, $\sigma_e=\sigma_{\rm SDSS}$)',fontsize=20)
plt.title(r'S21 prior (N = 445, $\sigma_{\rm e}=\sigma_{\rm OGLE}$)',fontsize=20)
plt.xlim(1,1000)
plt.grid()
plt.ylim(0.1,1.6)
plt.legend(fontsize=13,loc=4,ncol=2)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.show()

'''
for i in range(len(rho)):
    print(np.round(rho[i],1),
          np.round(np.mean(tau2,axis=1)[i,3]/np.mean(tau,axis=1)[i,3]-1,3))

